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Abstract 

Gurrently, most of the basic mechanisms governing tumor-immune system interactions, in combination 
with modulations of tumor-associated vasculature, are far from being completely understood. Here, we 
propose a mathematical model of vascularized tumor growth, where the main novelty is the modeling 
of the interplay between functional tumor vasculature and effector cell recruitment dynamics. Parame¬ 
ters are calibrated on the basis of different in vivo immunocompromised Ragl“/“ and wild-type (WT) 
BALB/c murine tumor growth experiments. The model analysis supports that tumor vasculature nor¬ 
malization can be a plausible and effective strategy to treat cancer when combined with appropriate 
immuno-stimulations. We find that improved levels of functional tumor vasculature, potentially medi¬ 
ated by normalization or stress alleviation strategies, can provide beneficial outcomes in terms of tumor 
burden reduction and growth control. Normalization of tumor blood vessels opens a therapeutic win¬ 
dow of opportunity to augment the antitumor immune responses, as well as to reduce the intratumoral 
immunosuppression and induced-hypoxia due to vascular abnormalities. The potential success of normal¬ 
izing tumor-associated vasculature closely depends on the effector cell recruitment dynamics and tumor 
sizes. Furthermore, an arbitrary increase of initial effector cell concentration does not necessarily imply 
a better tumor control. We evidence the existence of an optimal concentration range of effector cells for 
tumor shrinkage. Based on these findings, we suggest a theory-driven therapeutic proposal that optimally 
combines immuno- and vaso-modulatory interventions. 

Keywords: tumor-immune system interactions, murine tumor growth experiments, vascular normaliza¬ 
tion, immuno-modulatory interventions, mathematical modeling. 



2 


Major Findings 

We propose a tumor-effector cell recruitment model, where the main novelty is the low dimensional mod¬ 
eling of the complex interplay between functional tumor-associated vasculature and effector cell recruit¬ 
ment dynamics. Improved tumor vascularization, either via normalization or due to a microenvironmental 
stress alleviation strategy, is predicted to open a therapeutic window of opportunity in the treatment of 
cancer. Normalizing tumor vasculature as a potential therapeutic target closely depends on the immune 
recruitment dynamics and tumor sizes. Moreover, arbitrary low or high initial concentrations of effector 
cells can be detrimental to clinical outcomes, which evidences the existence of an optimal concentration 
range of effector cells that is crucial to obtain tumor control. 


Introduction 


Suppression of tumor angiogenesis has been recognized for more than four decades as a therapeutic 
target to treat solid tumors, as well as a complementary method of cancer prevention [l]-[^. Tradi¬ 
tionally, antiangiogenesis strategies attempt to inhibit new blood vessels formation within the tumor 
microenvironment, while at the same time reducing as much as possible the functionality of existing 
tumor-associated vasculature. A wide array of drugs to inhibit tumor angiogenesis, such as bevacizumab 
(Avastin, Genentech/Roche), a ligand-trapping monoclonal antibody against VEGF, and the inhibitors 
of multiple protein kinase targets, sorafenib (Nexavar, Bayer) and sunitinib (Sutent, Pfizer), has been 
successfully used in the clinic to treat solid tumors |^. However, different preclinical and clinical trials 
reveal that angiogenesis inhibitors alone have limited efficacy and only offer a modest survival benefit 
in a reduced cancer patient population [^[^. The clinical benefits of antiangiogenesis treatments are 
not only limited in terms of tumor shrinkage, vasculature destruction and patient survival, but also are 
restricted by transient effects, insufficient efficacy or development of treatment resistance [^|^[^. There 
are increasing evidences that the antitumor activity of most antiangiogenic drugs only became clinically 
significant in combination with conventional therapeutic modalities such as radiotherapy, chemotherapy 
immunotherapy [£ 11 


or 


Tumor-associated vasculature exhibits opposing effects on tumor growth, invasion and progression 
[^|^[^. On the one hand, blood vessels are required for delivery of oxygen and nutrients to the tumor, 
which promotes its growth, and they provide malignant cells with a way to spread, invade and metastasize 
into other healthy areas of the body (7|H2- 14 . On the other hand, blood vessels reoxygenate the tumor. 


which increases its radiosensitivity, and improve delivery of cytotoxic drugs to the tumor bulk as well as 
infiltration of immune cells into the tumor parenchyma [^[^[^[^. In view of these opposing effects, 
the right treatment of a cancer patient relies on a mechanistic understanding and quantitative analysis of 
the cancer-mediated processes involved. This demands for a mathematical model of tumor growth that 
considers the interplay between functional tumor vasculature and immune cell recruitment dynamics. 

There are several mathematical models of tumor growth, where some forms of the immune system 

Different immuno-oncology studies reveal that tumors can survive 

as inferred from clinical data 


dynamics were included 17- 21 


in a microscopic and undetectable dormant state 19 20 


22 


equilibrium can be disrupted by sudden events affecting the immune system, 
develops several strategies to circumvent the action of the immune system 


23 . This 


24 


25 


Indeed, the neoplasm 
which could result in 


tumor evasion and regrowth 25 . Sustained oscillations by the immune system have been observed both in 


its healthy state and pathological situations ^26[27 . An elegant mathematical model of the dynamics of 


immunogenic tumors revealed that tumor-immune system interactions can produce cyclic behaviors 28 


In fact, the presence of an immune component in mathematical models has been described to be crucial 
for reproducing observed phenomena such as tumor dormancy and oscillatory growth 22 28 -30 . This 


was complemented by a detailed discussion of the importance of including an immune component in 
tumor growth models, see 29 and references therein, and was subject to several reviews 17 31-35 . 
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Some of the models mainly focus on spatial aspects described either by partial differential equations or 
by cellular automatas or by realistic mechanical cell interactions [Mp3p6p7l 


while others focus on non- 


spatial dynamics described by ordinary differential equations |38| . Mathematical modeling provides 


a valuable theoretical framework to perform in silico experiments, as well as to verify assumptions and 
suggests modifications of existing theories. Despite years of research, and the vast amount of reported 
mathematical models, there are still many unanswered questions regarding the critical components of 
the tumor microenvironment and their role in the associated immuno-oncological mechanisms. However, 
models were rather successful in optimizing therapeutical protocols based on quantitative incorporation 
of experimental or patient data 29 39 45 . The work presented here, follows this latter approach. 

We intend to generate novel mechanistic insights not only on the interactions between the immune 
system and growing tumors, but also on the role of the functional tumor-associated vasculature. The 
main goal is to use the better understanding to suggest optimized treatment strategies. For that purpose, 
we combined a mathematical model of radially symmetric tumor growth 46 49 with an immune cell 
model 28 to propose a tumor-effector cell recruitment model. The main novelty is that the impact of 
functional tumor vasculature on the effector cell recruitment dynamics is included. Previous models of 
tumor-immune system interactions did not directly incorporate the effect of tumor vascularization, and 
thus underestimated the relevance of this crucial mechanism for tumor growth and treatment outcomes. 
Based on well-supported biological assumptions, we derive that tumor vasculature normalization opens 
a therapeutic window of opportunity for tumor growth control. The potential benefit of a vascular 
normalization strategy relies on a specific immune recruitment dynamic into the tumor microenvironment. 
We also find that an unlimited increase on the initial effector cell concentration does not necessarily 
imply a better tumor control. The proposed model is capable of providing a more complete view of 
vascularization effects on the immune responses to tumor growth, and it is therefore suited for analysis 
and prediction of treatment strategies. 


Quick Guide to Model Equations and Assumptions 


46 


A radial tumor growth model 

We consider a mathematical model of radial tumor growth which was originally proposed in 
also 1^ 49 . The tumor is assumed as an incompressible fluid flowing through a porous medium, where 
tissue elasticity is simplified. The tumor-host interface is considered to be sharp, and cell-to-cell adhesive 
forces are modelled as a surface tension at that interface. The tumor expands as a mass whose growth is 
governed by a balance between cell birth (mitosis) and death (apoptosis and necrosis). The mitotic rate 
within the tumor is assumed to be linearly-dependent on the nutrient concentration (oxygen, glucose, 
etc.) and is characterized by its maximal value Am at the tumor-host interface. The death rate 
is considered uniform within the tumor and assumed constant in time. The concentration of nutrients 
obeys a reaction-diffusion equation in the tumor volume, where nutrients are supplied from the functional 
tumor vasculature and consumed by the tumor cells at a uniform consumption rate (see the derivation 
of the mathematical model provided in the Supplementary Material for further details). 

The resulting dimensional mathematical model for the tumor radius R{t) dynamics, under the as¬ 
sumption of radial symmetry, is given by 


f = 1(A„B - A,)ii+ A„(I - B)L, 

This model given by Eq. 0 can be interpreted as the temporal evolution of the average tumor radius, 
since radial symmetric growth is not a realistic behaviour. Notice that, multiplication of Eq. 0byi?2 
describes the time evolution of the tumor volume (see the Supplementary Material for more details). The 
non-negative and dimensionless parameter B represents the net effect of tumor-associated vasculature on 
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the tumor radius dynamic. Moreover, Ld is an intrinsic scale representing the average length of nutrient 
gradient, i.e. supply, diffusion and consumption. 


A tumor-effector cell recruitment model considering the impact of functional 
tumor-associated vasculature 


The proposed tumor-effector cell recruitment model is a combination of the radial tumor growth model 
given by Eq. 0 and an effector cell model initially introduced in [^, see also (^|^. In addition, we 


consider the impact of functional tumor vasculature on the tumor growth and immune response dynamics. 
The system variables are the average tumor radius R(t) and effector cell concentration E(t) in the tumor 
vicinity. The model is formulated as a system of ordinary differential equations (ODEs) given by 


"I = ^(AmB - A,)fl + A„(l - B)Ln - cERm.o). (2) 

f]E d 3 

= rj^^E-{d,R^f{R,a)+do)E + <7, (3) 

where the time coordinate t on the system variables has been omitted for notational simplicity. The non- 
dimensionalized version of the model is provided in the Supplementary Material. In words, the Eqs. 0 
and 0 can be respectively read as follows 

{change of tumor radius} = {net vascular tumor growth} + {avascular tumor growth} - {death of tumor 
cells due to effectors}. 


{change of effector cell concentration} = {immune recruitment} - {inactivation + death of effector cells} 
+ {background rate of immune effector recruitment}. 


The expression f{R,a) = 


G [0,1], for (a G [0,1), is a phenomenological scaling function that 


models the infiltration of effector cells into the tumor bulk through the existing functional tumor vascu¬ 
lature. This function modulates the terms related to tumor-effector cell interactions, such as killing of 
tumor cells due to effectors and inactivation of effectors due to their antitumor activity. As a particular 
case, we assume that a = where 0 < B < 1 represents the net effect of functional levels of the 

tumor-associated vasculature. In the limit of avascular tumor growth 5 = 0, such tumor-effector cell 
interactions take place only at the tumor surface. At the other extreme, for 5 = 1 effector cells can 
potentially interact with any cancer cell within the tumor bulk. Notice that such scaling functions have 
been also considered in the classical von Bertalanffy approach and more recently in allometric models 


51 


Alternatively, we can view the encountering of effector and tumor cells mediated by fractal vasculature 
as a percolation process within the tumor mass. In this context, the term f{R,a) scales the average 
effector-tumor cell encountering rate in a vascularized tumor. In particular, the scaling term should be 
proportional to the corresponding percolation cluster mass of system size R and fractal dimension 3(a, 


where the latter depends on the structure of tumor-associated vasculature 52 


Model assumptions 

1. The temporal evolution of the average tumor radius is considered, where invasive and diffusive 
tumor mechanisms are not taken into account. 

2. The death rate of tumor cells reflects the lump effect of apoptotic and necrotic processes. 

3. Innate immunity or base immune surveillance is represented as a minimum presence of active effector 
cells at any time, even in the absence of cancer cells. 
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4. Effector cells are recruited depending on the concentration of tumor cells according to Michaelis- 
Menten dynamics. 


5. The efficacy of immune killing depends on the ability of effector cells to penetrate the tumor bulk via 
the functional tumor-associated vasculature 


10 [n]. With improved vascularization, the effectors 


kill tumor cells not only on the surface of the tumor, but also further inside. 


6. Effector cells die with a constant rate and get inactivated in dependence on their antitumor activity. 


Materials and Methods 

Experimental data of tumor growth in mice 

We now describe the murine experiments of in vivo tumor growth considered. BALB/c mice were pur¬ 
chased (Janvier, Erance). Ragl“/“ mice were purchased from The Jackson Laboratory (C.129S7(B6)- 
Ragl^^^^®^/J). lga~^~ mice were obtained from Michael Reth, Ereiburg, Germany (Cd79a^^^^^^). 
All recombinant mice were bred at the HZI and all experiments were performed with female, 8-12 weeks 
old mice if not stated differently and under approval of LAVES (Niedersachsisches Landesamt fiir Ver- 
braucherschutz und Lebensmittelsicherheit); Permission: 33.9-42502-04-12/0173. Unlike wild-type (WT) 
BALB/c mice, immunocompromised Ragl“/“ mice cannot produce functional T- and B-cells, i.e. have 
no adaptive immune responses. 

CT26 (ATCC CRL-2638) cells were cultured in IMDM supplemented with 10% ECS, 100 U/ml peni¬ 
cillin and 100 yug/ml streptomycin, 50 /iM 2-mercaptoethanol, 2 jaM L-glutamine and maintained at 37 °C 
and 5% C02- The El All (H-2^) cell line is a murine fibrosarcoma that expresses /3-galactosidase (/3-gal) 
and was obtained by transduction of spontaneously transformed BALB/c fibroblast cell line El with the 
LBSN retroviral vector [^. Tumors CT26 and El All were set by injecting 5 x 10^ cells in 100 /il PBS 
subcutaneously (s.c.). Growth was monitored by caliper and the mean tumor radius evolution between 
days 4 to 9 was observed. Volume was calculated as U = 4/37r(hu;^)/8, where h = height and w = width. 


Estimation of model parameters 

The experimental data corresponding to immunocompromised Ragl“/“ and wild-type (WT) BALB/c 
mice are analyzed in two different phases. In the first phase (days 0 to 4) exponential avascular tumor 
growth is assumed 54 . Using the tumor radius at those early times of growth, the net proliferation 
rate of cancer cells is estimated. The tumor radius evolutions in the second phase of growth (days 4 
to 9) are then used to determine the remaining model parameter values. Because of uncertainties due 
to the finite number of observations, consisting in 15 and 20 individual experimental trials respectively, 
we have implemented a bootstrapping procedure for each dataset considered. The goal was to obtain 
reliable estimates of the summary statistics, where a resampling procedure has been performed. A set 
of bootstrap samples is artificially generated and represents the natural variations in the experimental 
datasets 55 . More precisely, we have considered non-parametric bootstrap samples, i.e. samples that 


have been generated from the original experimental datasets by drawing with replacements. The default 
size managed, also called cloud size, has been 100 000 bootstrap samples per dataset, where the resulting 
mean and standard deviation of such samples have been considered for parameter calibration. Another 
crucial factor to obtain precise model fitting is to consider different initial parameter values. Gonsequently, 
we started the fitting procedure from a large number of different random initial conditions, as well as 
by randomly perturbing the model parameter set to be estimated. Then, the solution with the lowest 
residual variance was selected, which thus provides the best fit. 
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Early exponential tumor growth phase (days 0 to 4) 

For a precise understanding of the early tumor growth phase predicted by our mathematical model, 
we work under the small tumor radius assumption and Taylor expand the r.h.s. of Eq. 0 accordingly. 
Keeping terms up to the first nonlinear term, we obtain the following approximation 




(4) 




The very early tumor growth phase, dictated by the first term of this approximation, is always of 
exponential nature and does not depend on the parameter 5, i.e. the tumor vascularization effects, as 
reported also in [4^49|54| . Therefore, this initial phase depends exclusively on the net proliferation rate of 
cancer cells, Xp = (Am — A^)- Considering the initial tumor radius in the immunocompromised Ragl“/“ 
mice experiments of 0.40 mm, which corresponds to a tumor spheroid of 5 x 10^ cells of about 10 /im in 
diameter, and the tumor radius at day 4 of 2.0 mm, see Fig.m A), we find that Xp ^ 1.20T0.1 day“^. This 
value is in agreement with previous estimates of growth rate (early stages) for murine CT26 spheroids 56 


Later, the initial exponential tumor growth gets affected by tumor-associated vascularization, which starts 
playing a significant role when the two terms of the r.h.s. of Eq. 0 have similar orders of magnitude. 
This occurs around a critical tumor radius Rb derived from the condition 


-(Am — Xa)Rb 


Am(1 B) 
451/1) 




(5) 


Since (Am — A^) = Ap > 0 is necessary for the initial tumor growth, Eq. 0 requires 0 < 5 < 1, 
where B = 0 represents an avascular tumor and B ^ 1 a fully vascularized tumor. This yields 


Rb 


Ld 


15(Am — Xa) 
Am(1 - B) 


( 6 ) 


Eq. 0 provides an estimate of the tumor size at which nonlinearities, involving the vascularization 
mechanisms, start significantly influencing the initial exponential tumor growth. The effect of the non¬ 
linear terms is to saturate the growth, i.e. dR/dt = 0, when the tumor radius reaches Rb- Since only 
exponential growth is observed for small tumor radii, angiogenesis should be activated before the critical 
tumor radius is reached and suggests that Rb is an upper bound for the initiation of angiogenic processes. 
The value of Rb can be predicted provided the values of the parameters in Eq. 0 are known. Assuming 
time-invariant values of the characteristic mitotic Am and death A^ rates of cancer cells, and the relation 
Xp ^ 1.20 ± 0.1 day“^ together with the physiological plausible value Am = 1/18 h ^ 1.34 day“^ for the 
mouse colon carcinoma cell line (or CT26 murine tumor cells) (5^[57| , yields A^ ~ 0.14±0.1 day“^. The 
parameters B and Lb^ i-c. the functional tumor vasculature and intrinsic scale representing the average 
length of nutrient gradient, remain to be determined, which are estimated from the growth phase ranging 
between days 4 to 9. 


Linear tumor growth phase (days 4 to 9) 

From the experimental data corresponding to immunocompromised Ragl“/“ and (WT) BALB/c mice 
at the stages of growth ranging between days 4 to 9, we observe an approximately linear evolution of the 
tumor radius, see Fig. This observation is consistent with the long-term behavior of growth, where 
linear radial evolution is prominent as reported in [^. In Eqs. 0-0 this assumption is equivalent to 
assigning B = Xa/Xm- In particular, linear radial tumor growth in immunocompromised Ragl“/“ mice 
requires c = 0 in Eq. (|^, because Ragl“/“ mice have no adaptive anticancer immune responses. Thus, 
the first and last terms on the r.h.s. of Eq. Q vanish and allows to determine Lb ^ 0.29±0.02 mm. This 
value is consistent with the well-known characteristic nutrient diffusion length that ranges between 0.2 mm 
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and 0.3 mm [^[^ . In Fig. [^A) the linear fitting results for the Ragl“/“ mice dataset are represented. 
The estimated parameter values allow for the calculation of Rb ~ 1.12 ±0.08 mm > Lb from Eq. ^ for 
a tumor in the early phase of growth which begins to be influenced by the tumor-associated vasculature. 

2Qp8 |5Q|[6Q| [^ 


The following parameter values are taken from 
K = 2.72 mm^ and a = 0.13 x 10^ cells day“^ 


do = 0.37 day di = 0.01 mm day 


-1 


CT26 growth in Ragl-/- mice CT26 growth in WT BALB/c mice 



Time (days) Time (days) 


Figure 1: Linear fitting of the tumor radius evolution from murine experimental data. (A) 

Tumor growth in immunocompromised Ragl“^“ and (B) wild-type (WT) BALB/c mice at different time 
points are shown together with the linear fitting results of the model. The ranges of experimental data 
and standard deviations of the bootstrap samples are also shown. 


Radial tumor growth can be also assumed linear for the (WT) BALB/c mice dataset, see Fig. Bb), 
which infers dR/dt = Vl ~ 0.29 ±0.05 mm/day constant in the time interval between days 4 to 9, and in 
agreement with growth velocity measured previously for murine CT26 tumors 
Eq. § becomes 


62 . Then, ioi R^ Lb ^ 


^{XmB - \a)R + Am{1 - B)Ld - cER^ = Vl, (7) 

where (^ tanh(A/Li 3 ) “ ^ for 0 < 5 < 1 it follows that since <C 1. 

Considering that R{t) = Vl(^ ~ ^o) ± Eq. 0 provides the following evolution of the effector cell 
concentration E{t) required for linear radial tumor growth 


E{t) 


—Vl ± ^{XmB — XA)R{t) ± Am(1 — B)Lb 
cR{t)B 


( 8 ) 


The above formula allows for the exact quantification of the initial concentration of effector cells Eq, 
as well as for the reduction of the number of free model parameters that enter in the fitting procedure, 
since the net proliferation Xp and death rates of cancer cells, as well as Lb have been estimated from 
the experimental data of immunocompromised Ragl“/“ mice tumor growth. Accordingly, only three free 
model parameters are left in Eqs. Q-(|^, i.e. functional tumor vasculature 5, immune cell recruitment 
rate r and death rate of tumor cells due to effectors c, which are estimated from the (WT) BALB/c mice 
experiments imposing the expression ^ during parameter calibration. Considering variations of the 
initial conditions, and randomly perturbing the parameter set to be estimated, the resulting parameter 
values are provided in Tab. In Fig. [^B) the linear fitting results for the experimental data of tumor 
growth in (WT) BALB/c mice are represented. 











Parameters Values 
c 0.03 ± 0.002 cells-i 

^ n /I ^ n ni 


Sources 




Xa 

Ld 

B 

Eo 


day 


0.29 ± 0.02 mm 

0.18 ± 0.01 

3.28 ± 1.07 X 10^ cells 


estimated 

estimated 


K 

2.72 mm^ 

20 

28 

50 

61 

di 

0.01 mm“^ day“^ 

^0 

28 

50 

61 

do 

0.37 day“^ 

20 

28 

60 

61 

(J 

0.13 X 10^ cells day“^ 

20 

28 

50 

61 

Am 

1.34 day“^ 

56 

57 




estimated 

estimated 

estimated 

estimated 


Table 1: Model parameters. Parameter values not found in the literature were estimated from in 
vivo experimental data of tumor growth in immunocompromised Ragl“^“ and wild-type (WT) BALB/c 
mice. 


Results 


The dual effect of functional tumor-associated vasculature 

The proposed mathematical model given by Eqs. ©-<© is used to investigate the effect of functional 
vasculature on the gross tumor growth. Tumor-associated vascularization enhances the supply of oxygen 
and nutrients, and therefore favours tumor growth and progression 12 14 . At the same time blood 


vessels facilitate the infiltration of effector cells into the tumor bulk, which potentially exterminate cancer 
cells 15 16 . These opposing effects of tumor vascularization suggest the existence of an optimal tumor 


growth control, such that a quantitative modulation of functional vasculature can be considered as a 
suitable therapeutic target against cancer. 

In order to gain insights into the impact of functional vasculature on the short-term tumor evolution, 
we vary the model parameter B in Eq. (§. The radial growth rate dR/dt^ which indicates whether the 
tumor radius grows (positive), decays (negative) or remains stable (zero), is monitored. Starting from 
different initial concentrations of effector cells Eq, Eig.[^ shows the rate of tumor growth in dependence 
on the functional levels of vascularization B and initial tumor sizes Rq. Above a specific threshold 
concentration of effector cells, a degree of vascularization exists that maximizes tumor growth. Moreover, 
above this threshold increasing vascular functionality reduces the tumor growth rate, since it allows a 
high number of effector cells to penetrate into the tumor parenchyma, see diagrams for Eq > 20 x 10^ cells 
in Eig.[^ Tumor shrinkage dR/dt < 0 is possible for very low functional vascularization and large tumor 
radii, irrespective of the number of effector cells in the tumor vicinity. However, for small tumors and 
low initial concentrations of effector cells, reduction in the tumor burden cannot be obtained for deficient 
vascularization. Moreover, there exist critical Eq values at which tumor shrinkage becomes also possible 
for well-vascularized tumors irrespective of their sizes, see diagram for Eq = 30 x 10^ cells in Eig. 
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Effector Cells Effector Cells Effector Cells 

(Eo = 3.3 X 105 cells) (Eo= 10 x 10^ cells) (Eo= 15 x 10^ cells) 



Effector Cells Effector Cells Effector Cells 

(Eo= 20 X 105 cells) (Eo = 25 x 105 cells) (Eo = 30 x 105 cells) 



Figure 2: Short-term growth rate dependence on tumor size Rq and functional vascularization 

B. The radial tumor growth rate (dR/dt) is plotted against initial tumor radius Rq and functional levels 
of tumor-associated vasculature B for different initial concentrations of effector cells Eq. Remaining 
parameters in Eq. Q are as in Tab. 

A more systematic investigation of the dependence of long-term model dynamics given by Eqs. (H)-® 
on the set of initial conditions (Eq, Rq) is conducted through a tumor-immune phase portrait analysis, 
see Eig. Eor functional vascular levels of B = 0.50, long-term tumor control is possible, in the sense 
of damped oscillations, for all values of Rq considered, see Eig. |^A). Tumor dormancy is described by 
decaying oscillations converging to a steady-state, that once obtained is maintained for a long period of 
time. However, at larger values of Rq the amplitude of the oscillations in E{t) and R{t) increases and 
reaches uncontrolled growth above a threshold for Rq which determines the tumor fate. Eor a higher 
tumor-associated vascularization, e.g. B = 0.70, long-term tumor control cannot be achieved for any set 
of values of the model parameters considered, see Eig. [^B). In contrast, for well-vascularized tumors, 
e.g. B = 0.95, a window of opportunity for tumor control is predicted, see Eig. [^C). As far as Rq stays 
below a threshold value, tumor evolution can be controlled. At high and low vascularization levels, tumor 
eradication is not possible for sufficiently large tumors, revealing the relation between tumor sizes and the 
potential therapeutic benefit of targeting the tumor vascular network. The bifurcation analysis on model 
parameters B and r, i.e. the functional tumor vasculature and immune recruitment rate, is provided in 
the Supplementary Material for further details. 
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Figure 3: Tumor-immune phase portrait diagrams and long-term characteristic dynamics. 

(Left) Phase portrait diagrams from top to bottom representing the tumor radius R and concentration 
of effector cells E time evolutions for different sets of initial conditions (colored curves) and levels of 
functional tumor vascularization B. The nullclines, i.e. zero-growth isoclines, of the dynamical system 
(black lines) are also represented. Each arrow provides information about the short-term (local) evolution 
of the corresponding initial conditions. (Right) Panels showing the corresponding time evolution of R 
(red lines) and E (blue lines) for Eq = 20 x 10^ cells and different values of the initial tumor radius Rq 
for each value of B considered in the corresponding phase portrait diagrams on the left. The immune 
recruitment rate was r = 0.59 day“^ and remaining parameters are as in Tab. [l] 

Critical radius and control range of effector cell concentration determine tu¬ 
mor fate 

Fig. [^reveals information about the existence of different dynamical regimes of the tumor-effector cell 
recruitment system described by Eqs. (©-(§ for the set of parameter values considered. Based on the 
analysis of the critical fixed points and their positions in the phase space (see the phase space and critical 
fixed points analysis provided in the Supplementary Material), we obtain that a spiral and a saddle point 
are present. At the spiral point the tumor radius oscillates around a small tumor size, i.e. a dormancy- 
associated equilibrium point. This implies a concentration range of effector cells that keep tumor under 
control for a period of time that depends on the point stability. Moreover, there exists a critical tumor 
radius Rs^ which coincides with the saddle point, that for tumors of equal or larger sizes always grow 
uncontrollably and the immune system is completely suppressed. However, below this threshold tumors 
are always controlled, or at least for some period of time in a transient dormant state, provided the 
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immune responses are in a specific range. The questions that arise are to (i) analytically evaluate the 
critical tumor size Rs that for Rq > Rs implies uncontrollable growth and (ii) estimate the range of 
initial effector cell concentration Eq that allows for long-term tumor control. 

Critical radius for tumor dormancy 

These results evidence that the saddle point position in the phase space separates uncontrolled tumor 
growth and immune-induced dormant states, as well as delineates the potential and limits of therapeutic 
benefits of immuno- and vaso-modulatory interventions. For large well-vascularized tumors, we derive an 
analytical expression that approximates the critical tumor radius Rs corresponding to an upper bound 
estimation of the saddle point (see the estimation of the critical radius for tumor dormancy provided in 
the Supplementary Material for further details), which is given by 


Rs = 


Sac 


di{XM — Xa) 


2(do - r) 
di 


(9) 


We demonstrate that for values of the tumor radius equal or higher than Rs^ the eigenvalues corre¬ 
sponding to the radial time evolution are positive. This implies that for radii beyond Rs, the tumor grows 
uncontrollably evading immune responses. The nonlinear system of equations §-(§ was numerically 
solved using the Trust-Region Dogleg method 63 . Fig.j^A) shows the dependence of Rs on the immune 
recruitment rate r obtained by means of model simulations. Rs represents an upper bound of the critical 
tumor size that determines the long-term tumor fate, and as r increases the simulated and estimated 
values get closer. For R < Rs long-term or at least transient tumor dormancy can be obtained, while 
R > Rs results in tumor escape from the immune system surveillance. In the situation of tumors with 
low functional vascular levels, Rs can be directly estimated by means of model simulations. 


A 



(D 

DC 


B AE (x 105 cells) 



Vasculature (B) 


Figure 4: Dependence of the estimated critical tumor radius Rs^ the tumor radius at the 
unstable saddle point and control range of effector cell concentration AE on the immune 
recruitment rate r. (A) The tumor radius at the unstable saddle point in dependence on the parameter 
r derived from model simulations (blue dashed line) and compared to the approximation of Rs given by 
Eq. ([]) (red solid line) for a degree of functional tumor vascularization 5 = 0.95. (B) AE with respect to 
model parameters r and B. Remaining parameters are as in Tab. 


Control range of effector cell concentration regulates tumor evolution 

The efficacy of infiltrating effector cells to control the tumor growth is derived from the properties of the 
spiral point, which is always obtained for tumors with size lower than the critical threshold Rs. Tumor 
dormancy, and then neoplasia control, takes place when the tumor radius stays close to the spiral point 
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for a long period of time. Therefore, the attraction range of the spiral point with respect to the E-axis 
provides useful information about the number of effector cells required for long-term tumor control. 

We claim that for arbitrary high or low initial concentrations of effector cells Eq, immune-induced 
tumor dormancy is not guaranteed. To illustrate the counterintuitive loss of tumor control for high or low 
initial numbers of effector cells, we focus at the phase portrait diagram of Fig. [^B). There we observe 
that when the amplitude of the E{t) oscillations exceeds a certain range, the tumor escapes from the 
immune system surveillance. Mathematically this is translated to the existence of a concentration range 
of effector cells AE that is crucial for the control of tumors with size smaller than the critical radius Rs- 
Details of AE calculation are provided in the Supplementary Material. 

Fig. [^B) shows the dependence of the control range of effector cell concentration AE on functional 
tumor-associated vasculature B and immune recruitment rate r. In particular, we observe that AE 
increases for high enough immune recruitment rates and increasing tumor vascularization. However, at 
low values of r < 0.54 day“^, AE becomes smaller as B increases. The latter is intuitive considering that 
low immune recruitment rates could not provide enough effector cells to exploit the enhanced potential 
infiltration offered by increasing functional tumor vasculature. Therefore, the pro-tumoral effect of the 
functional vascular network overcomes its antitumoral action. 


Therapeutic potential of immuno- and vaso-modulatory interventions 

Immune-mediated tumor dormancy depends not only on the effectiveness of the immune system responses, 
but also on the functional tumor-associated vasculature. In consequence, we investigate the impact 
of model parameters B and r, i.e. the functional tumor vasculature and immune recruitment rate, 
on the overall therapeutic potential for long-term (5-years) model predictions. Therapeutic potential 
TP{p) : ftp [0,1] is interpreted as the average tumor control, i.e. tumor radius stays bounded, with 
respect to a certain regime of model parameter p G {H, r} given by 


TP{p) = 


J H[Rs{p) - R{p,t*)]dp 
J H[R{p,t*)]dp 

Op 


( 10 ) 


where t* = 1825 days (5-years) is the target time, Rs{p) is the tumor radius at the saddle point and i^(-) 
is a Heaviside step function. The effect of the initial tumor radius Rq and effector cell concentration Eq 
have been also investigated. It should be noted that a 5-years tumor-free patient is typically considered 
in a complete remission stage, i.e. permanent absence of disease, and thus cured. 

Any increase in the initial tumor radius Rq always reduces the therapeutic potential of the immune 
system to control tumor growth, see Fig. However, an optimal concentration of effector cells for 
long-term tumor control is predicted irrespective of the degree of functional tumor vascularization B. 
This observation implies that an arbitrary increase in the initial concentration of effector cells does not 
necessarily result in a greater therapeutic benefit. It is counterintuitive that high values of Eq do not 
enhance long-term tumor control. A large initial number of effector cells can induce a great initial tumor 
burden reduction that makes the tumor escape immune control by an undercritical immune stimulus. 
In mathematical terms, the tumor radius time evolution crosses the separatrix, escapes from the limit 
cycle and enters into the unlimited growth dynamic regime (see for instance the green curve in phase 
portrait diagram of Fig. [^B), as well as the phase space and critical fixed points analysis provided in the 
Supplementary Material). 
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0.55 < Vasculature (B) < 0.85 C 0.85 < Vasculature (B) <1.0 D 




— Tumor Radius Ro = 2.5 mm m Tumor Radius Ro = 3.0 mm 
Tumor Radius Ro = 3.5 mm — Tumor Radius Ro = 4.0 mm 


Figure 5: Therapeutic potential with respect to the initial concentration of effector cells Eq 
and functional tumor vascularization B for different values of the initial tumor radius Rq 
and immune recruitment rate r. The parameter r varies between 0.50 day“^ and 0.70 day“^ for 
each value of Eq considered. The y-axis represents the ratio of 5-years tumor control situations over all 
possible cases computed for the values of r and B in the ranges considered. Remaining parameters are 
as in Tab. [H 

In addition to the initial concentration of effector cells, the immune recruitment dynamics play an 
important role in long-term tumor control, see Fig. Large values of the immune recruitment rate r 
enhance the tumor control, irrespective of the tumor vascularization regime B. However, for high levels 
of functional tumor vasculature, i.e. H > 0.55, large tumors with Rq > 4 mm are almost uncontrollable 
at any immune recruitment rate r, see Fig. [^C,D). We also demonstrate that the initial concentration 
of effector cells required for long-term tumor control mostly depends on functional tumor vascularization 
(see the phase space and critical fixed points analysis provided in the Supplementary Material). 
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0.50 0.55 0.60 0.65 0.70 
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— Tumor Radius Ro = 2.5 mm — Tumor Radius Ro = 3.0 mm 
A Tumor Radius Ro = 3.5 mm # Tumor Radius Ro = 4.0 mm 


Figure 6: Therapeutic potential with respect to the immune recruitment rate r and functional 
tumor vascularization B for different values of the initial tumor radius Rq and concentration 
of effector cells Eq. The parameter Eq varies between 10 x 10^ cells and 30 x 10^ cells for each value 
of r considered. The y-axis represents the ratio of 5-years tumor control situations over all possible cases 
computed for the values of Eq and B in the ranges considered. Remaining parameters are as in Tab. 


A theory-driven therapeutic proposal 


Model predictions provide new insights that can be used to suggest and design novel therapeutic strate¬ 
gies against tumor growth. The aim is to propose alternative or complementary ideas for less toxic cancer 
treatments compared to conventional therapeutic modalities, such as surgery, radiotherapy or chemother¬ 
apy. The latter cancer therapies are related to immunosuppressive effects that disallows the intention to 
exploit immune system responses. 

An appropriate tumor profiling via medical imaging techniques and tissue biopsy is essential. Imaging 
allows for the estimation of the initial tumor size Rq and its comparison with the critical radius Rs 
that potentially determines the long-term tumor fate. For calculation of the precise Rs value, the net 
proliferation rate of tumor cells could be extrapolated from biopsy samples by following the methodology 
proposed in 64 . Moreover, from the same pathology data the initial concentration of effector cells Eq 
can be estimated, e.g. by counting the number of CD8+ anti-tumor effector T cells. Then, the range of 
effector recruitment rate could be assumed comparable to the one estimated from our experimental data, 
see Tab. For Rq < Rs mid- or long-term controlled tumor growth is possible by clinically influencing 
the interplay between immune recruitment dynamics and functional tumor-associated vasculature, as 
suggested by the therapeutic potential results in Figs. and In the situation of Rq > Rs^ tumor 
control could not be obtained by the findings obtained from the proposed model, and complementary 
treatment modalities are initially required. 

Positive therapeutic responses, dR/dt < 0, can be achieved by increasing or decreasing the functional 
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tumor vasculature. The latter is possible via antiangiogenic drugs, such as bevacizumab, sorafenib or 
sunitinib. However, low levels of tumor vascularization trigger hypoxic-induced phenomena, and this 



opportunity in terms of long-term tumor control for different functional levels of tumor-associated vas¬ 
culature H, initial number of effector cells Eq and tumor sizes Rq. A negative correlation between initial 
tumor sizes and mid- or long-term control is predicted. For sufficiently small tumors, optimal Eq values 
range with higher therapeutic potentials exist. An increase of the immune recruitment rate results in a 
higher therapeutic potential and the window of opportunity for high functional levels of tumor vasculature 
still persists, see Fig.|^ Therefore, in addition to tumor vasculature normalization, applying a monitored 
immunotherapy protocol could further increase the expected therapeutic benefits. Fig. describes the 
complete theory-driven therapeutic proposal. 
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Figure 7: Theory-driven therapeutic proposal. Medical imaging and biopsy allow for estimating 
the initial tumor radius Rq and concentration of effector cells Eq^ as well as the immune recruitment 
r and tumor characteristic mitotic Am rates. If the tumor is sufficiently small, R < Rs^ a vasculature 
normalization strategy can be implemented to substantially improve tumor control. The potential benefit 
of normalizing the tumor vasculature is conducted through an increase of r and optimizing Eq. 


Conclusions 

In this work, we have explored the possibilities offered by mathematical modeling and computer simu¬ 
lations to gain insights into the interplay between functional tumor-associated vasculature and immune 
recruitment dynamics. We have combined a radially symmetric tumor growth and an effector cell model 
including the impact of vascularization on the immune system efficacy in killing tumor cells. The model 
parameter calibration was based on murine experiments of in vivo tumor growth. Model analysis revealed 
the existence of a tumor control window of opportunity when high functional levels of tumor vascula¬ 
ture are present. This can be biologically perceived since high vascularization allows for the penetration 
of effector cells into the tumor bulk and diminishes the adverse effects of hypoxia occurred due to low 
oxygenation. Improving tumor vascularization can be mediated by a normalization of blood vessel func¬ 
tionality [^[7|[^ [^[6^ . Such an approach is successful only if the tumor size is below a critical threshold 
that depends on tumor cell metabolic needs. Furthermore, we predict that an arbitrary increase in the 
initial effector cell concentration does not necessarily imply a direct enhancement of tumor control, thus 
optimal values are needed. According to our model predictions when the initial reduction of tumor size 
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exceeds a specific threshold, the tumor escapes immune control and relapses faster than the associated 
tumor-induced immune responses. Taking into account the dependency of tumor control on these various 
mechanisms, we have proposed an alternative less toxic therapeutic approach based on a personalized 
tumor profiling. This therapeutic proposal, based on normalizing the tumor vasculature, can be also 
potentially combined with chemotherapy and radiotherapy to enhance the therapeutic outcomes mm 


We conclude by pointing out a number of related future research directions, as well as the limitations 
of this work. The current results are based on murine experiments and therefore a parameterization on 
human data would be essential. Introducing tumor vascularization dynamics should be not only a natural, 
but also an insightful extension of the proposed mathematical model. Moreover, we need to further refine 
the modeling of the term /(i?, a) and elaborate the dependence of the exponent a on the functional tumor 
vasculature B. Additionally, we could extend the notion of parameter a to reflect the dynamic feedback 
of effector cell activity on local tumor density fluctuations as a time-evolving fractal dimension. At this 
state, vascularization B is considered as a static parameter in time. Making tumor vasculature dynamic 
would make sense to model further hypoxic effects, such as necrosis 73 or hypoxic-induced invasion 48 


Moreover, this will make possible to further investigate the dual effects of functional blood vessels on 
tumor growth dynamics when cytotoxic drugs are administered as part of chemotherapy. On the other 
hand, the immune system is much more complex than its description in the current model, involving a 
wide range of immune cell types and intricate interactions. In particular, the immune system is often 
regarded as acting to suppress tumor growth, however it is now clear that it can be both stimulatory and 
inhibitory 74 . Recent advances have indicated that tumor-associated macrophages actively promote all 
aspects of tumor initiation, progression and development 75 . Therefore, modeling the interplay between 


macrophages and tumor cells, while taking into account the functional tumor-associated vasculature and 
effector cell recruitment dynamics, may highlight new targets to develop novel anticancer strategies. 
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Supplementary Material 

Derivation of the mathematical model of tumor growth 

We consider a non-necrotic tumor which volume growth results from a balance between cell mitosis and 
death, driven by the presence of nutrients, e.g. oxygen or glucose. In the absence of inhibitor chemical 
species, the spatio-temporal dynamic of the nutrient concentration cr(x, t) is modelled by the quasi-steady 
react ion-diffusion equation given by 


o = L)vV + r, 

where D is the diffusion coefficient and F is the rate at which nutrients are supplied to the tumor volume 
Q{t). The quasi-steady assumption is well supported by the observation that the diffusion time scale for 
oxygen or glucose is much lower 1 minute) than the tumor cells doubling time 1 day). The rate F 
incorporates all sources and sinks in the tumor volume, and is based on the following phenomenological 
assumptions: 

(1) Nutrients are homogeneously supplied by the vasculature at a rate F^ = —Xb {cf — cfb)^ where is 
the uniform nutrient distribution in the blood and is uniform. 

(2) Nutrients are consumed by tumor cells at a rate Xa. This yields a rate F given by 

F = -Xb {cr - (Jb) - Xa . 

The tumor is modelled as an incompressible fluid which velocity field u in satisfies the continuity 
equation given by 


V•u = Ap, 

where Ap is the net proliferation rate of tumor cells that leads to volume growth (or shrinkage). This 
formulation rests on the additional assumptions described below: 

(3) The tumor is modelled as a unique homogeneously distributed phenotype, meaning that all tumor 
cells behave in exactly the same manner. 

(4) Tumor expansion exclusively depends on the net cell proliferation, and invasive processes, e.g. cell 
migration or diffusion, are not explicitly included. 

(5) The model assumes that the density of tumor cells is constant and homogeneous within the tumor bulk. 

(6) The net cell proliferation rate is assumed as follows 

Ap = - Xa , 

g-oo 

where is the nutrient concentration outside the tumor volume, and the mitotic Am and death rates 
A^ of tumor cells are uniform. 

(7) The velocity is assumed to obey the Darcy’s law, i.e. porous media flow, given by 


u = —/iVP, 
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where /i is a (constant) cell motility parameter and P(x, t) is the pressure inside the tumor that is assumed 
to satisfy the Laplace-Young boundary condition at the interface, which corresponds to the following as¬ 
sumption: 

( 8 ) Cell-to-cell adhesive forces are modelled by a surface tension 7 at the tumor boundary. By introducing 
the intrinsic length scale Ld = j {Xb + we obtain an intrinsic relaxation time scale = 

We consider these length and time scales to non-dimensionalize the mathematical model 
that can be rewritten, using bar-notation for dimensionless quantities, in terms of the modified nutrient 
concentration a and pressure p defined by 


and 


P 


where the parameter 


7 

Ld 



^Am( 1 - P)(l - cr) + {Xa - XmB) 


X • xA 

) 


B 

~ (7^ XbPX 

represents the net effect of vascularization. 

Then, by using algebraic manipulations, the original dimensional problem can be reformulated in 
terms of two non-dimensional decoupled problems as follows 

V2d-d = 0, 

= 1 ; 

and 


= 0 , 


/_N 1 /^ ^ (X • X)S 

(p)e = /^ - “ AmP)—^—, 

in a d-dimensional tumor separated from the host tissue by the interface Y of local curvature k that 
evolves with the normal velocity = n • (u)e 5 and n being the outward normal to Y. 

When considering evolution of a three-dimensional tumor that remains radially symmetric, problems 
above have analytical solutions that lead to the evolution equation for the dimensionless tumor radius R 
given by 


dR _ 1 ^ 1 

~ 

In particular, we have considered the dimensional version of the previous equation for the tumor 
radius R = LdR evolving with respect to time t = I/Xr. 

Tumor volume evolution 

It is instructive to compare the proposed mathematical model of tumor growth with other modelling 
approaches. In particular, we compare it with the well-known logistic growth model which has been 
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extensively used in previous studies. To do that, we multiply both sides of Eq. (1) in the manuscript by 
3i?^, and for large radii, i.e. R ^ Ljj, we obtain the following equation 

dTi 

3i?2— = {XmB - Xa)R^ + 3Am(1 - B)LdR^ 

where we have that 


lim 

R^oo 


( _ ^ _ 

\tdinh{R/LD) 


Ld \ 
R J 


= 1 . 


Now, assuming that the tumor mass is M = pV oc R^ and p is the tumor density per volume element, 
then we obtain an equation for the tumor mass given by 


— = {XmB - Xa)M + 3Am(1 - B)LdM^/\ 

where the term denotes the mass on the surface of volume V. In particular, for B = 0, i.e. an 

avascular tumor, the above equation becomes similar to the well-known von Bertallanfy equation 


— = 3XmLdM^/^ - XaM, 

which represents the growth of a body that is supplied with nutrients only from its surface. The last 
equation is a kind of generalized logistic growth. In order to derive the classical logistic growth equation 
we need to restrict ourselves to a two-dimensional system. As a matter of fact, the main problems with 
the logistic model are: 


1. It is valid only for very large tumors, but it fails to capture the dynamics of small tumors. 

2. It is true only for fixed nutrient consumption within the tumor volume, i.e. no nutrient diffusion. 

3. It is valid only for two-dimensional tumor evolution. 


Finally under the same assumption of R Ljj and for a fixed concentration of effector cells Eq, we 
have that 


— = {XmB - Xa)M + 3Xm(,1 - B)LdM^/^ - cEqM^ . 

The last term represents the “depth” of tumor-effector cell interactions with respect to the tumor- 
associated vasculature. Obviously, for the avascular limit B = 0 such term indicated only surface inter¬ 
actions, while for B = 1 effectors can kill cancer cells throughout the tumor bulk. 


Model non-dimensionalization 


The process of non-dimensionalization allows for (i) scaling the system of variables according to their 
characteristic quantities, (ii) writing the system by reducing as much as possible the number of model 
parameters, and (hi) identifying their appropriate units. To do that, we assume as reference length size 
the model parameter and the initial number of effector cells as reference concentration Eref^ for the 
variables R and E respectively. The diffusion length Ld ranges in the interval [0.2 mm, 0.3 mm] 5^59[76 , 
and the effector cell characteristic concentration E^ef is at the order of magnitude 10^ cells. The latter 
estimate is justified since the characteristic length scale of the system is at the order of 1 mm, and given 
that cells are commonly assumed with a diameter between 10 pm and 20 pm 57 , then for a volume of 
1 mm^ the concentration is at the order of magnitude 10^ cells. Moreover, the time scale is proportional 
to tumor cell deactivation, i.e. t' = {EEref)t^ where the parameter c' is provided below. 
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Thus, the non-dimensional system of ordinary differential equations for the new variables R' = R/Lu 
and E' = E jE^ef is given by 




^ - c'E' 


R 


/B 




R' 


rE'-{d[- 


R' 


dE' 

IF ~ ^ K'+R^^ y 


+ djQ j E' -|- 


where the non-dimensional model parameters are summarized as follows 




c' Ere f ^ 


A'4 = 77 


Aa 


c' Ere f ^ 


c' Eref ’ 


K' = 


K 


d[ 


diL 


(2 + B) 


c'E, 


ref 


d'o 


dp 


c' Eref ’ 


„ • 

ref 


It should be noted that the inactivation rate of effectors by encounters with tumor cells interestingly 
depends on the parameter 5, i.e. the functional tumor vasculature, by means of the dimensionless 
parameter d'l. 


Phase space and critical fixed points 

The nullclines, i.e. zero-growth isoclines, as well as critical fixed points corresponding to different model 
parameter values of the tumor-effector cell recruitment dynamical system {R-E axis), are shown in Fig.[^ 
The intersections of the nullclines, i.e. the curves along which dR/dt = 0 and dE/dt = 0, represent the 
critical fixed points of the system. Fig. shows that a spiral and a saddle point are present, which are 
represented by S and U on each phase diagram respectively. At the spiral point the time evolution of R 
oscillates around a small tumor size, and the separatrix defines sharp boundary between tumor growth 
and shrinkage depending on the initial conditions, see also Fig. 3 in the manuscript. There exists a critical 
tumor size that coincides with the saddle point and above which tumors always grow uncontrollably and 
the immune system is completely suppressed. Below this threshold, tumors are always controlled when 
the spiral point is stable, or at least for some period of time if it is unstable, provided the immune 
responses are in a specific range. Therefore, the mathematical model is capable of explaining both tumor 
dormancy and escape from immunosurveillance. 
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r = 0.54 days-1 / B = 0.50 
S = (1.98, 20.34) / U = (3.37, 21.05) 



r = 0.54 days-1 / B = 0.95 
S = (2.14, 25.26) / U = (2.83, 25.75) 



r = 0.64 days-1 / B = 0.50 r = 0.64 days-i / B = 0.95 

S = (1.60, 20.31) / U = (4.23, 21.69) S = (1.62, 25.29) / U = (3.65, 25.48) 



Figure 8: Nullclines and critical fixed points. The nullclines dE/dt = 0 (blue lines) and dR/dt = 0 
(red lines) are plotted for different values of model parameters r and i.e. immune recruitment rate 
and functional tumor-associated vasculature respectively. The intersections of the nullclines represent 
the spiral and saddle points of the system, denoted by S and U respectively. Remaining parameters are 
as in Tab. 1 of the manuscript. 

A quick glance at Fig. [^reveals that the distance between the critical fixed points depends on the model 
parameter values. Fig. [^represents the position of the critical fixed points on the R-E axis for different 
values of functional tumor vascularization B and immune recruitment rate r. The variations of the tumor 
radius corresponding to the spiral point are always smaller, ranging between 1.59 mm to 2.08 mm, than 
for the unstable saddle point, varying from 2.82 mm to 5.80 mm, see Fig.|^A,B) respectively. For B ^ 1 
and small values of r, the spiral and saddle points approach each other in the R-axis. Accordingly, this 
means that as the immune recruitment rate decreases and functional tumor vascularization increases, the 
tumor becomes less controllable, i.e. the region for favorable initial conditions in terms of long-term tumor 
control is lower. Conversely, at high immune recruitment rates and low tumor vascular functionality, the 
tumor can be controlled in a range of initial conditions, which is consistent with expectations. The initial 
concentration of effector cells E only weakly depends on the immune recruitment rate r, but mostly on 
the level of functional tumor vascularization 5, see Fig. me ,D). Concentrations of effector cells associated 
with the spiral and the saddle points both increase with increasing tumor vasculature, and approximately 
to the same extent, which corresponds to a parallel shift on the F^-axis of both critical fixed points, see 
Fig. me, D). 
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Figure 9: Position of the critical fixed points in the phase space. Position dependence of the 
spiral and saddle points on different values of model parameters r and i.e. immune recruitment rate 
and functional tumor-associated vasculature respectively. (A,B) Tumor radius and (C,D) effector cell 
concentration corresponding to the spiral (A,C) and saddle points (B,D). Remaining parameters are as 
in Tab. 1 of the manuscript. 


Bifurcation analysis 

In order to explore the critical parameter values at which the qualitative long-term dynamics of tumor 
evolution changes, we generated one-parameter bifurcation diagrams for the functional tumor vasculature 
B and immune recruitment r, with the tumor radius plotted as a function of such parameters. Stable 
and unstable solution branches are represented by solid and dashed lines respectively. The Hopf and 
limit bifurcation points are indicated by H and LP respectively. A Hopf point is a local bifurcation 
in which a fixed point of a dynamic system loses stability as a pair of complex conjugate eigenvalues 
of the linearization around the fixed point cross the imaginary axis of the complex plane, i.e. when the 
eigenvalues are a complex conjugate pair whose real part changes sign. The bifurcation analysis has been 
conducted by a prediction-correction continuation algorithm using the software MATCONT 
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For functional levels of the tumor-associated vasculature in a range around H = 0.75 and bound by the 
Hopf bifurcation points Hi and 772, the spiral fixed point is a repellor, see Fig.[^ Consequently, for B in 
this range and for the different values of r, uncontrolled tumor growth occurs. Out of this critical range, 
i.e. for lower and higher values of H, the spiral fixed point is an attractor, where long-term tumor control 
is reached. Not only low levels of functional vasculature allows for long-term tumor control, but there 
exists a therapeutic window of opportunity for well-vascularized tumors, see Fig. With increasing 
immune recruitment rate r, the distance between the Hopf bifurcation points decreases, implying that 
the observed window of opportunity increases, while at the same time the global region of unfavourable 
situations decreases. 
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Figure 10: Bifurcation diagrams with respect to the functional tumor-associated vascu¬ 
lature B, One-parameter bifurcation diagrams for different values of the immune recruitment rate 
r = {0.54, 0.57, 0.60, 0.64} day“^ in reading order. The labels Hi and H 2 represent Hopf bifurcation 
points, while the solid and dashed lines describe the stable and unstable solution branches respectively. 
The concentration of effector cells E corresponding to Hi and H 2 are also provided. Remaining param¬ 
eters are as in Tab. 1 of the manuscript. 


In the case of one-parameter bifurcation diagrams with respect to the immune recruitment rate r, 
we observe a limit and a Hopf bifurcation point. When tumor vascular functionality B 


see Fig. 11 


increases, the stable branch is gradually reduced and the Hopf bifurcation point moves from low to high 
values of r. However, there exists a value of H, at which this behaviour is inverted and the stable branch 
gets larger again. The range of r values associated with the stable branch for B = 0.95 is similar to 
the one obtained for low values of H, for instance B = 0.40, see Fig. [m In other words, even for a 
tumor with a highly functional vascular network, long-term control can be achieved for a wide range of 
immune recruitment rates. This finding suggests an alternative therapeutic strategy to the conventional 
antiangiogenic treatments. 
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Figure 11: Bifurcation diagrams with respect to the immune recruitment rate r. One- 
parameter bifurcation diagrams for different values of functional tumor-associated vasculature B = 
{0.40, 0.60, 0.80, 0.95} in reading order. The labels H and LP represent a Hopf and a limit bifur¬ 
cation point respectively. The solid and dashed lines describe the stable and unstable solution branches 
respectively. The concentration of effector cells E corresponding to H and LP are also provided. Re¬ 
maining parameters are as in Tab. 1 of the manuscript. 

Estimation of the critical radius for tumor dormancy 

We derive an analytical expression that approximates the R-E axis position of the saddle point, which 
corresponds to an upper bound estimation. Ld and B ^ 1 imply in Eq. (2) in the manuscript that 
^ ^ 0 and tanh (R/Ld) 1, as well as that • These simplifications and the quasi-steady 

state approximation of Eq. (2) yields 

— ^a)R — 

which implies that 


^ _ 2(Am — Aa) 

“ 3^ ■ 

Arguing similarly for Eq. (3) in the manuscript and under the same assumptions, we have that 
1 and where 

f R^ \ 

rE — Idi ——h do j -E + (j = 0, 

and then we obtain that 


a 


E = - 


r - {di^ + do) 















26 


Combination of the previous equations for E results in an approximation of the value of R corre¬ 
sponding to the critical tumor radius, i?-axis position of the saddle point, given by 


Rs 


3ac 

di{XM — Xa) 


2{do - r) 
di 


Estimation of control range of effector cell concentration 

We define the control range of effector cell concentration AE as the maximum difference of the ^-axis 
spiral point position and the number of effectors £’* given by the separatrix that defines sharp boundary 
between tumor dormant states and uncontrollable growth. Combinations of estimated AE with the 
position of the critical fixed points in the phase space for a given set of model parameters provide 
information about the ability of the immune responses to induce long-term tumor control. Values of 
initial concentration of effector cells Eq out of AE will not result in long-term tumor dormancy, which 
suggests limitations of therapeutic immuno-modulatory interventions. 


References 

[1] J. Folkman, “Anti-angiogenesis: new concept for therapy of solid tumors,” Ann Surg, vol. 175, 
pp. 409-416, Mar. 1972. 

[2] R. Jain, “Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combi¬ 
nation therapy,” Nat Med^ vol. 7, pp. 987-989, Sept. 2001. 

[3] Z. Miao, J. Feng, and J. Ding, “Newly discovered angiogenesis inhibitors and their mechanisms of 
action,” Acta Pharmacol Sin, vol. 33, pp. 1103-1111, Sept. 2012. 

[4] A. Argyriou, E. Giannopoulou, and H. Kalofonos, “Angiogenesis and anti-angiogenic molecularly 
targeted therapies in malignant gliomas,” Oncology, vol. 77, pp. 1-11, May 2009. 

[5] J. Ebos and R. Kerbel, “Antiangiogenic therapy: impact on invasion, disease progression, and metas¬ 
tasis,” Nat Rev Clin Oncol, vol. 8, pp. 210-221, Mar. 2011. 

[6] G. Jayson, D. Hicklin, and L. Ellis, “Antiangiogenic therapy-evolving view based on clinical trial 
results,” Nat Rev Clin Oncol, vol. 9, pp. 297-303, Eeb. 2012. 

[7] R. Jain, “Normalization of tumor vasculature: An emerging concept in antiangiogenic therapy,” 
Science, vol. 307, pp. 58-62, Jan. 2005. 

[8] R. Jain, “Normalizing tumor microenvironment to treat cancer: bench to bedside to biomarkers,” J 
Clin Oncol, vol. 31, pp. 2205-2218, June 2013. 

[9] D. Duda, R. Jain, and G. Willett, “Antiangiogenics: the potential role of integrating this novel 
treatment modality with chemoradiation for solid cancers,” J Clin Oncol, vol. 25, pp. 4033-4042, 
Sept. 2007. 

[10] Y. Huang, J. Yuan, E. Righi, W. Kamoun, M. Ancukiewicz, J. Nezivar, M. Santosuosso, J. Mar¬ 
tin, M. Martin, E. Vianello, P. Leblanc, L. Munn, P. Huang, D. Duda, D. Eukumura, R. Jain, and 
M. Poznansky, “Vascular normalizing doses of antiangiogenic treatment reprogram the immunosup¬ 
pressive tumor microenvironment and enhance immunotherapy,” Proc Natl Acad Sci USA, vol. 109, 
pp. 17561-17566, Oct. 2012. 



27 


[11] Y. Huang, S. Goel, D. Duda, D. Fukumura, and R. Jain, “Vascular normalization as an emerging 
strategy to enhance cancer immunotherapy,” Cancer Res, vol. 73, pp. 2943-2948, May 2013. 

[12] N. Ferrara and R. Kerbel, “Angiogenesis as a therapeutic target,” Nature, vol. 438, pp. 967-974, 
Dec. 2005. 

[13] S. Weis and D. Cheresh, “Tumor angiogenesis: molecular pathways and therapeutic targets,” Nat 
Med, vol. 17, pp. 1359-1370, Nov. 2011. 

[14] R. Farnsworth, M. Lackmann, M. Achen, and S. Stacker, “Vascular remodeling in cancer,” Oncogene, 
vol. 33, pp. 3496-3505, Aug. 2014. 

[15] W. Fridman, F. Pages, C. Sautes-Fridman, and J. Galon, “The immune contexture in human tu¬ 
mours: impact on clinical outcome,” Nat Rev Cancer, vol. 12, pp. 298-306, Mar. 2012. 

[16] M. Junttila and F. de Sauvage, “Influence of tumour micro-environment heterogeneity on therapeutic 
response,” Nature, vol. 501, pp. 346-354, Sept. 2013. 

[17] N. Bellomo and M. Delitala, “From the mathematical kinetic, and stochastic game theory to mod¬ 
elling mutations, onset, progression and immune competition of cancer cells,” Phys Life Rev, vol. 5, 
pp. 183-206, Dec. 2008. 

[18] J. Arciero, T. Jackson, and D. Kirschner, “A mathematical model of tumor-immune evasion and 
sirna treatment,” Discrete Continuous Dyn Syst Ser B, vol. 4, pp. 39-58, Feb. 2004. 

[19] A. Matzavinos, M. Ghaplain, and V. Kuznetsov, “Mathematical modelling of the spatio-temporal 
response of cytotoxic t-lymphocytes to a solid tumour,” Math Med Biol, vol. 21, pp. 1-34, Mar. 2004. 

[20] A. d’Onofrio, “A general framework for modeling tumor-immune system competition and im¬ 
munotherapy: Mathematical analysis and biomedical inferences,” Physica D: Nonlinear Phenomena, 
vol. 208, pp. 220-235, Sept. 2005. 

[21] M. Robertson-Tessi, A. El-Kareh, and A. Goriely, “A mathematical model of tumor-immune inter¬ 
actions,” J Theor Biol, vol. 294, pp. 56-73, Feb. 2012. 

[22] G. Garavagna, A. d’Onofrio, P. Milazzo, and R. Barbuti, “Tumour suppression by immune system 
through stochastic oscillations,” J Theor Biol, vol. 265, pp. 336-345, Aug. 2010. 

[23] G. Koebel, W. Vermi, J. Swann, N. Zerafa, S. Rodig, L. Old, M. Smyth, and R. Schreiber, “Adaptive 
immunity maintains occult cancer in an equilibrium state,” Nature, vol. 450, pp. 903-907, 2007. 

[24] D. Pardoll, “Does the immune system see tumors as foreign or self?,” Annu Rev Immunol, vol. 21, 
pp. 807-39, 2003. 

[25] G. Dunn, L. Old, and R. Schreiber, “The three es of cancer immunoediting,” Annu Rev Immunol, 
vol. 22, pp. 329-60, 2004. 

[26] J. Stark, G. Ghan, and A. George, “Oscillations in the immune system,” Immunol Rev, vol. 216, 
pp. 213-231, Apr. 2007. 

[27] B. Goventry, M. Ashdown, M. Quinn, S. Markovic, S. Yatomi-Glarke, and A. Robinson, “Grp iden¬ 
tifies homeostatic immune oscillations in cancer patients: a potential treatment targeting tool?,” J 
Transl Med, vol. 7, p. 102, 2009. 



28 


[28] V. Kuznetsov, 1. Makalkin, M. Taylor, and A. Perelson, “Nonlinear dynamics of immunogenic tumors: 
Parameter estimation and global bifurcation analysis,” Bull Math Biol^ vol. 56, pp. 295-321, Mar. 
1994. 

[29] L. de Pillis, W. Gu, and A. Radunskaya, “Mixed immunotherapy and chemotherapy of tumors: 
modeling, applications and biological interpretations,” J Theor Biol^ vol. 238, pp. 841-862, Feb. 
2006. 

[30] A. d’Onofrio, F. Gatti, P. Gerrai, and L. Freschi, “Delay-induced oscillatory dynamics of tumourim- 
mune system interaction,” Math Comput Models vol. 51, pp. 572-591, Mar. 2010. 

[31] R. Araujo and D. McElwain, “A history of the study of solid tumour growth: the contribution of 
mathematical modelling,” Bull Math Biol, vol. 66, pp. 1039-1091, Sept. 2004. 

[32] H. Byrne, T. Alarcon, M. Owen, S. Webb, and P. Maini, “Modelling aspects of cancer dynamics: a 
review,” Philos Trans A Math Phys Eng Sci, vol. 364, pp. 1563-1578, June 2006. 

[33] T. Roose, S. Ghapman, and P. Maini, “Mathematical models of avascular tumor growth,” SIAM 
Rev, vol. 49, pp. 179-208, May 2007. 

[34] R. Eftimie, J. Bramson, and D. Earn, “Interactions between the immune system and cancer: a brief 
review of non-spatial mathematical models,” Bull Math Biol, vol. 73, pp. 2-32, Jan. 2011. 

[35] K. Wilkie, “A review of mathematical models of cancer-immune interactions in the context of tumor 
dormancy,” Adv Exp Med Biol, vol. 734, pp. 201-234, Aug. 2013. 

[36] G. Schaller and M. Meyer-Hermann, “Multicellular tumor spheroid in an off-lattice voronoi/delaunay 
cell model,” Phys Rev E, vol. 71, pp. 051910-1-16, 2005. 

[37] G. Schaller and M. Meyer-Hermann, “A modelling approach towards an epidermal homoeostasis 
control,” J Theor Biol, vol. 247, pp. 554-573, Aug. 2007. 

[38] R. Sachs, L. Hlatky, and P. Hahnfeldt, “Simple ode models of tumor growth and anti-angiogenic or 
radiation treatment,” Math Comput Model, vol. 33, pp. 1297-1305, June 2001. 

[39] H. Byrne, “Dissecting cancer through mathematics: from the cell to the animal model,” Nat Rev 
Caneer, vol. 10, pp. 221-230, Mar. 2010. 

[40] H. V. Jain and M. Meyer-Hermann, “The molecular basis of synergism between carboplatin and 
abt-737 therapy targeting ovarian carcinomas,” Caneer Res, vol. 71, pp. 705-715, 2011. 

[41] D. Gorwin, G. Holdsworth, R. Rockne, A. Trister, M. Mrugala, J. Rockhill, R. Stewart, M. Phillips, 
and K. Swanson, “Toward patient-specific, biologically optimized radiation therapy plans for the 
treatment of glioblastoma,” PLoS ONE, vol. 8, p. e79115, Nov. 2013. 

[42] H. Kempf, H. Hatzikirou, M. Bleicher, and M. Meyer-Hermann, “In silico analysis of cell cycle 
synchronisation effects in radiotherapy of tumour spheroids,” PLoS Comput Biol, vol. 9, p. el003295, 
Nov. 2013. 

[43] J. Alfonso, N. Jagiella, L. Nunez, M. Herrero, and D. Drasdo, “Estimating dose painting effects in 
radiotherapy: a mathematical model,” PLoS ONE, vol. 9, p. e89380, Eeb. 2014. 

[44] J. Alfonso, G. Buttazzo, B. Garcfa-Archilla, M. Herrero, and L. Niihez, “Selecting radiotherapy dose 
distributions by means of constrained optimization problems,” Bull Math Biol, vol. 76, pp. 1017- 
1044, May 2014. 



29 


[45] H. V. Jain, A. Richardson, M. Meyer-Hermann, and H. M. Byrne, “Exploiting the synergy between 
carboplatin and abt-737 in the treatment of ovarian carcinomas,” PLoS One, vol. 9, p. e81582, 2014. 

[46] H. Greenspan, “On the growth and stability of cell cultures and solid tumors,” J Theor Biol, vol. 56, 
pp. 229-242, Jan. 1976. 

[47] H. Byrne and M. Chaplain, “Growth of nonnecrotic tumours in the presence and absence of in¬ 
hibitors,” Math Biosci, vol. 130, pp. 151-181, Dec. 1995. 

[48] V. Cristini, J. Lowengrub, and Q. Nie, “Nonlinear simulation of tumor growth,” J Math Biol, vol. 46, 
pp. 191-224, Mar. 2003. 

[49] H. Hatzikirou, A. Chauviere, J. Lowengrub, J. De Groot, and V. Cristini, “Effect of vascularization 
on glioma tumor growth,” in Modeling Tumor Vaseulature (T. L. Jackson, ed.), pp. 237-259, Springer 
New York, 2012. 

[50] L. de Pillis, A. Radunskaya, and C. Wiseman, “A validated mathematical model of cell-mediated 
immune response to tumor growth,” Caneer Res, vol. 65, pp. 7950-8, Sept. 2005. 

[51] A. Herman, V. Savage, and G. West, “A quantitative theory of solid tumor growth, metabolic rate 
and vascularization,” PLoS ONE, vol. 6, no. 9, p. e22973, 2011. 

[52] D. Stauffer and A. Aharony, Introduetion To Pereolation Theory. Taylor & Erancis, 1994. 

[53] P. Paglia, C. Chiodoni, M. Rodolfo, and M. Colombo, “Murine dendritic cells loaded in vitro with 
soluble protein prime cytotoxic t lymphocytes against tumor antigen in vivo.,” J Exp Med, vol. 183, 
no. 1, pp. 317-322, 1996. 

[54] A. Bru, S. Albertos, J. Garcfa-Asenjo, and I. Brii, “Pinning of tumoral growth by enhancement of 
the immune response,” Phys Rev Lett, vol. 92, pp. 1-4, June 2004. 

[55] A. Davison and D. Hinkley, Bootstrap Methods and their Applieation. Cambridge: Cambridge Uni¬ 
versity Press, 1997. 

[56] K. Alessandri, B. Sarangi, V. Gurchenkov, B. Sinha, T. KieBling, L. Eetler, E. Rico, S. Scheming, 
C. Lamaze, A. Simon, S. Geraldo, D. Vignjevic, H. Domejean, L. Rolland, A. Eunfak, J. Bibette, 
N. Bremond, and P. Nassoy, “Cellular capsules as a tool for multicellular spheroid production and 
for investigating the mechanics of tumor progression in vitro,” Proe Natl Aead Sei USA, vol. 110, 
pp. 14843-14848, Aug. 2013. 

[57] M. Delarue, E. Montel, O. Caen, J. Elgeti, J. Siaugue, D. Vignjevic, J. Prost, J. Joanny, and 
G. Cappello, “Mechanical control of cell flow in multicellular spheroids,” Phys Rev Lett, vol. 110, 
p. 138103, Mar. 2013. 

[58] H. Erieboes, X. Zheng, C. Sun, B. Tromberg, R. Gatenby, and V. Cristini, “An integrated computa¬ 
tional/experimental model of tumor invasion,” Caneer Res, vol. 66, pp. 1597-1604, Eeb. 2006. 

[59] V. Cristini, H. Erieboes, X. Li, J. Lowengrub, P. Macklin, S. Sanga, S. Wise, and X. Zheng, “Non¬ 
linear modeling and simulation of tumor growth,” in Seleeted topies in eaneer modeling: Genesis, 
evolution, immune eompetition, and therapy. Modelling and Simulation in Seienee, Engineering, and 
Teehnology (N. Bellomo, M. Chaplain, and E. de Angelis, eds.), ch. 6, pp. 113-82, Boston, MA USA: 
Birkhauser, 2008. 

[60] B. Su, W. Zhou, K. Dorman, and D. Jones, “Mathematical modelling of immune response in tissues,” 
Comput Math Methods Med, vol. 10, no. 1, pp. 9-38, 2009. 



30 


[61] A. d’Onofrio, U. Ledzewicz, and H. Schattler, “On the dynamics of tumor-immune system in¬ 
teractions and combined chemo-and immunotherapy,” in New Challenges for Caneer Systems 
Biomedieine, pp. 249-266, Springer, 2012. 

[62] M. Seshadri, J. Spernyak, P. Maiery, R. Cheney, R. Mazurchuk, and D. Bellnier, “Visualizing the 
acute effects of vascular-targeted therapy in vivo using intravital microscopy and magnetic reso¬ 
nance imaging: correlation with endothelial apoptosis, cytokine induction, and treatment outcome,” 
Neoplasia, vol. 9, pp. 128-35, Feb. 2007. 

[63] J. Nocedal and S. Wright in Numerieal Optimization, Springer Series in Operations Research and 
Financial Engineering, Springer-Verlag New York, 2006. 

[64] A. Hyun and P. Macklin, “Improved patient-specific calibration for agent-based cancer modeling,” 
J Theor Biol, no. 317, pp. 422-424, 2013. 

[65] A. Hoeben, B. Landuyt, M. Highley, H. Wildiers, A. Van Oosterom, and E. De Bruijn, “Vascular 
endothelial growth factor and angiogenesis,” Pharmaeologieal Reviews, vol. 56, pp. 549-580, Dec. 
2004. 

[66] D. Quail and J. Joyce, “Microenvironmental regulation of tumor progression and metastasis,” Nat 
Med, vol. 19, pp. 1423-1437, Nov. 2013. 

[67] S. Goel, D. Duda, L. Xu, L. Munn, Y. Boucher, D. Eukumura, and R. Jain, “Normalization of the 
vasculature for treatment of cancer and other diseases,” Physiol Rev, vol. 91, pp. 1071-1121, July 
2011. 

[68] T. Stylianopoulos, J. Martin, V. Chauhan, S. Jain, B. Diop-Erimpong, N. Bardeesy, B. Smith, 
C. Eerrone, E. Hornicek, Y. Boucher, L. Munn, and J. R. K., “Causes, consequences, and remedies 
for growth-induced solid stress in murine and human tumors,” Proe Natl Aead Sei USA, vol. 109, 
pp. 15101-15108, Aug. 2012. 

[69] E. Manning, J. Ullman, J. Leatherman, J. Asquith, T. Hansen, T. Armstrong, D. Hicklin, E. Jaf- 
fee, and L. Emens, “A vascular endothelial growth factor receptor-2 inhibitor enhances antitumor 
immunity through an immune-based mechanism,” Clin Caneer Res, vol. 13, pp. 3951-3959, July 
2007. 

[70] J. Hamzah, M. Jugold, E. Kiessling, P. Rigby, M. Manzur, H. Marti, T. Rabie, S. Kaden, H. Grone, 
G. Hammerling, B. Arnold, and R. Gauss, “Vascular normalization in rgs5-deficient tumours pro¬ 
motes immune destruction,” Nature, vol. 453, pp. 410-414, May 2008. 

[71] R. Shrimali, Z. Yu, M. Theoret, D. Chinnasamy, N. Restifo, and S. Rosenberg, “Antiangiogenic 
agents can increase lymphocyte infiltration into tumor and enhance the effectiveness of adoptive 
immunotherapy of cancer,” Caneer Res, vol. 70, pp. 6171-6180, Aug. 2010. 

[72] V. Agrawal, S. Maharjan, K. Kim, N. Kim, J. Son, K. Lee, H. Choi, S. Rho, S. Ahn, M. Won, 
S. Ha, G. Koh, Y. Kim, Y. Suh, and Y. Kwon, “Direct endothelial junction restoration results in 
significant tumor vascular normalization and metastasis inhibition in mice,” Oneotarqet, vol. 5, no. 9, 
pp. 2761-77, 2014. 

[73] P. Macklin and J. Lowengrub, “Nonlinear simulation of the effect of microenvironment on tumor 
growth,” J Theor Biol, vol. 245, no. 4, pp. 677-704, 2007. 

[74] K. de Visser, A. Eichten, and L. Coussens, “Paradoxical roles of the immune system during cancer 
development,” Nat Rev Caneer, vol. 6, pp. 24-37, Jan. 2006. 



31 


[75] N. Hao, M. Lii, Y. Fan, Y. Cao, Z. Zhang, and S. Yang, “Macrophages in tumor microenvironments 
and the progression of tumors,” Clin Dev Immunol^ vol. 2012, p. 11, June 2012. 

[76] H. Acker, Spheroids in eaneer researeh: methods and perspeetives. Springer-Verlag Berlin; New York, 
1984. 

[77] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “Matcont: A matlab package for numerical bifurca¬ 
tion analysis of odes,” ACM Trans Math Softw, vol. 29, pp. 141-164, June 2003. 



